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A variational approach is considered to calculate the free energy and the conformational 
properties of a polyelectrolyte chain in d dimensions. We consider in detail the case of 
pure Coulombic interactions between the monomers, when screening is not present, in 
order to compute the end-to-end distance and the asymptotic properties of the chain as 
a function of the polymer chain length N. We find R ~ A''(log A*')^ where ly — and 
A is the exponent which characterize the long-range interaction U cx The exponent 

7 is shown to be non-universal, depending on the strength of the Coulomb interaction. 
We check our findings, by a direct numerical minimization of the variational energy for 
chains of increasing size 2** < A*' < 2^^. The electrostatic blob picture, expected for small 
enough values of the interaction strength, is quantitatively described by the variational 
approach. We perform a Monte Carlo simulation for chains of length 2* < A*' < 2^^°. 
The non universal behavior of the exponent 7 previously derived within the variational 
method, is also confirmed by the simulation results. Non-universal behavior is found 
for a polyelectrolyte chain in d = 3 dimension. Particular attention is devoted to the 
homopolymer chain problem, when short range contact interactions are present. 



I. INTRODUCTION 

Charged polymers play a central role in theory and experiments. Polyelectrolytes are macromolecules 
containing ionizable groups. They dissociate in water to form charged groups and low-molecular- 
weight counter ions. It is important to remark that interactions between different polyions in a 
many chain system are of crucial importance for any realistic system with a low but finite overlap 
concentration. The problem of the single chain conformation is the one of the essential questions for 
polyelectrolyte solutions and has been under investigation in the recent years. This will be the topic of 
the present paper as well, but in a more general context. As common to most of the recent investiga- 
tion, we consider the limit of infinite dilution and zero salt concentration for polyelectrolytes. In this 
limit a polyelectrolyte can be represented as a connected sequence of charged and uncharged monomers 
in a "dielectric vacuum" that constitutes the solvent. ^From a scaling point of view, the properties of 
single polyelectrolyte chains, despite the long time debates are still not completely clear. Indeed, for the 
neutral polymer chain system, where a clear understanding of the scaling properties has been achieved 
0, the details of the local chain structure are not essential, while in the polyelectrolyte chain problem 
they could be important |^,^ . In particular the long range nature of the Coulomb interaction couples via 
the presence of the solvent, the short and long length scales of the system . Thus a scaling approach to 
polyelectrolytes results to be more involved than for neutral chains. The interaction between monomers, 
counterions and solvent molecules leads to the important phenomenon of counter ion condensation Q and 
a realistic description should take into account all the different characters of this scenario. On the other 
hand, a minimal or bare essential description of the system is required in order to achieve a quantitative 
description of the single polyelectrolyte chain problem. 

In this paper we will present a generalized form of des Cloizeaux Gaussian variational method |^,^ 
to study the scaling behavior and local properties of polyelectrolytes. Variational techniques have been 
shown to be well suited to describe a self- interacting polymer with long range repulsion Scaling 
properties of the end-to-end distance of the chain have been discussed for the long - ranged intra- molecular 
interaction at an arbitrary spatial dimensionality d. A recent numerical implementation of such a 
variational technique, for the case of Coulombic interaction (when A = d — 2), seems to be promising 
as well 1^. On the other hand, the validity of the variational approach for short range interacting 
homopolymer chains must be questioned |^-^ for natural and well known reasons. The most studied 
example in the context of polymer physics is the problem of the self avoiding chain. It was proven that the 
problem shows universality in the limit of infinitely long chains Q and can be treated by field theoretic 
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methods common to phase transitions. Indeed the renormahzation group (RG) theory manifests these 
statements and develops systematic methods to compute the scahng exponents, e.g. for the size of the 
chain ||^. The variational technique provides an alternative approach since most of the approximations 
can be controlled. Although, it is well known from elementary statistical mechanics that the variational 
technique cannot provide such powerful tool due to their mean - field character. Therefore this technique 
is expected to work in the mean field limit, i.e. at large dimensions (at the upper critical dimension) or 
for long range potentials. Thus it is well investigated and generally believed that the Gaussian variational 
approach is better suited for long-range interacting chains. 

Moreover, it has been suggested that in the study of random heteropolymer and directed 

polymer in a random media |12f| , a variational approach in replica space can be also successfully applied. 
We postpone these issues for future investigations. Here we will mainly concentrate on the pure long-range 
and short-range homopolymer problems within the Gaussian variational technique. 

In section II we introduce the Gaussian polymer model and calculate the variational free energy using 
a general Gaussian trial Hamiltonian in the spirit of des Cloizeaux and reference ||l^. We then 
critically review the asymptotic analysis given originally in reference |5|] and suggest a new one which 
leads to the Flory exponents corrected by the presence of logarithmic factors. In section III we describe 
a numerical technique in order to solve the Euler equations corresponding to the variational principle, 
by means of a new algorithm specially devised for this problem. An independent algorithm, that simply 
minimizes the variational free energy is shown to be quantitatively equivalent to the direct solution of the 
Euler equations. The numerical solution of the Euler equations is then carefully analyzed for polymer 
chains of increasing size. We show here that our new asymptotic analysis is in quantitative agreement 
with our numerical findings. We discuss connections with the short range homopolymer problem within 
the Gaussian Variational Principle and results for the end-to-end distance are obtained. In section 
IV Monte Carlo simulation results, obtained by means of a pivot algorithm, are presented, showing 
quantitative agreement with our variational analysis. The case of a self-interacting homopolymer chain, 
when an attractive two-body interaction between the monomers is considered together with a three body 
repulsive term, can also be studied within our variational approach. Preliminary calculations shows that 
the globular formation, expected from simple phenomenological arguments and due to the balancing 
between the two body and three body interactions, takes place. The instability of the globule state due 
to a Coulombic repulsion between the monomers is also described by our variational principle and will 
be presented elsewhere llj] . Conclusions and further prospects are summarized in section V. 



II. THE GAUSSIAN VARIATIONAL PRINCIPLE 



A. Basic equations 



First we start from a well known example for a Gaussian variational approach suitable for the investi- 
gations presented below. We introduce a simple problem to compute the free energy of a self-interacting 
polymer chain with Coulombic repulsion between the monomers, based on a discrete representation of a 
chain of N monomers. The method relies on the well known variational principle, where a Gaussian trial 
probability is proportional to the exponential of a quadratic form of the monomer coordinates rj- , 

Pv{rl,..,rN) = Zy^c^p{-Ho{rl,..,rN)}. (2.1) 
In eq.( ^.l| ) the trial Hamiltonian Hq is given by 

, N N 

Ho{rl,..,r-^) = -^^G^{^-J){n~r-)' (2.2) 



and Zy is a normalization constant determined by the condition 

Pv{rl,..,rN)d''rl...d''rN - 1. (2.3) 



Up to this point different forms of the Gaussian variational principle could be considered. We decided 
here to consider the variational form originally proposed by des Cloizeaux in the context of the self- 
avoiding homopolymer problem. The main essential features of the Gaussian var iational principle can 



be obtained within this choice of the variational form according to equation (2^). More sophisticated 
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forms of the variational kernel including a classical path fo and eventual anisotropy along parallel and 
perpendicular directions to the polymer backbone, can be also included [|l3|. By the appropriate choice 
of a periodic boundary conditions for the chain the quadratic form Hq can be diagonalized by introducing 
cyclic coordinates (Rouse modes) for the monomer positions 

N 

Pq ~ N~2 ^exp[ i2TTjq/N]'rj 

N 

= II exp[-i2TTjq/N]p'g. (2.4) 

q=l 

Indeed, due to cyclic invariance the Cartesian components of pq satisfy 

(pi'^P^f) = ls,,,Sqq,g],\27Tq/N), (2.5) 
where gJf^(2Trq/N) is positive and discontinuous function of k = 2'Kq/N and is related to the Gaussian 



propagator G of equation (2.2) via the following expression 



9N 



{2TTq/N) = ^ GN{n)[l ~ cos{2TiqnlN)\. (2.6) 



The variational free energy for a uniformly charged polymer of polymerization index N uiad dimensional 
space can be computed according to the usual Feynman inequality 

F<Fv = {H-Ho)o + Fo, (2.7) 



where Fq is the free energy of the Gaussian model defined by equation (2.2). In equation (2.7) H represents 
the full Hamiltonian of a charged polymer chain with a pure long - range form for the monomer-to- 
monomer interaction, 

Vi\r-l-r-}\) = \r-l-r-}\-\ (2.8) 
where X < d. The expression for the Hamiltonian reads 

, N N N 

i—l i—1 j—1 

where /3 — l^/A^ depends on the Bjerrum length l-g — {e'^ /ATreksT) and where the distance between 
charged monomers along the chain is given by A in units of the Kuhn segment b. A defines thus the 
charge fraction / — 1/A. 

It is straightforward to derive the following expression for the variational free energy Fy in the limit 
of infinitely long chain N ^ oo, 

Fy[gik)]/N = -Ll^ Y.^Kn)]-^ + 



\ 2 / n—1 



^ J^^ log g{k)dk - ^[1 + log(27r/d)], 



(2.10) 



where 



u r- ■^ 1 Z"^" (l-COs((^- j)A:))_,, ,^ ^ ^ , 



bN{^-J)^m-r;)^), (2.12) 



and 
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lim gN{kN) ^ gik), g{k) > 0. 



(2.13) 



The function g{k) satisfies the symmetry requirements 

gik)^gi-k), g{k + 27r)^g{k). (2.14) 
Minimization of the variational free energy Fy yields 

gik)^l^COs{k)-p(^^ r(^) 



2A2; r(f) 

^(l-cos(nfc))[6(n)]-^. (2.15) 



.2 

00 



X 

n=l 



Equations ( ^.1C| ) -(^.15| ) represent the generalization of the Gaussian variational principle for the long- 
range potential ( |2.8|) which was introduced first in refere nce [ zf. For the pure Coulombic case X = d—2, 
comparing the expression for the variational free energy ( ^.10 ) with the correspondent expression for the 
short range interaction problem |^ , 

, 00 
Fy[g{k)]/N^-b{l)+wY,[b{n)]-U 



2 

n=l 



^ J \ogg{k)dk^^[l + \og{2TT/d)], (2.16) 

leads to the following conclusion which holds within the Gaussian Variational Principle. All universal ex- 
ponents for the Coulombic case are related to the corresponding exponents for the short-range interaction 
case through the dimensional shift 

(2.17) 



The result of minimization of equation (2.16) reads 

00 

g{k) = 1 - cos(A:) - ^^(1 - cos(nfc))[6(n)]"^, (2.18) 

n=l 

which is of the same form like equation ( 2.15| ). 

Bef o re mo ving to the next section, where a careful numerical implementation of the Euler equations 
( 2.11 )-( p.l5| ) for polymer chains of increasing length TV is presented, we revise the asymptotic analysis 
of the Euler equations in the limit of small k 0, which was originally given by des Cloizeaux [^ . Th e 
behavior of g{k) at small k determines the asymptotic behavior of b{n) for n >> 1 via equation ( 2.11 ). 
Let us assume first that in this limit g{k) has the following form 

5(fc) ^ ^ < ^ < 1, (2.19) 

then the correlation function b{n), reads 

b{n) ~ bn^". (2.20) 

We will now proceed to analyze the Euler equations in the scaling limit of small k values within this 
simple ansatz. 



B. The asymptotic analysis - neutral chains 



To demonstrate several problems within the Gaussian v ariat ional principle we are going to rewise in 
this section the asymptotic analysis of the Euler equation (2.18) for neutral chains, which was originally 
given in the context of the self avoiding chain problem , and later extended to self interacting polymer 
chains with long-ranged Coulomb like interactions It has been suggested that due to the long ranged 
nature of the interactions between the monomers the Gaussian Variational Principle is particularly well 
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suited to describe the asymptotic behavior of polyelectrolyte chains. Special attention has also been 
devoted [p p^Jl5[ | to the case of screened Coulombic interaction, within the Debye Hiickel approximation. 
Let us consider the interaction term 



I{k) = - cos{nk))[b{n)Y 



(2.21) 



We can obtain the asymptotic expression for I{k) in the fc ^ limit, substituting the series with an 
integral Q and assuming that the dominant contribution would correspond to large values of n. The 
result can be written in the form 



/(fc) 



d + 2 
TTO 2 



2r((5)cos(f<5) 



(2.22) 



where 5 = v{d + 2). 

An alternative derivation can be obtained in the asymptotic regime n >> 1 evaluating exactly the 
series 



(2.23) 



at z = exp(iA:) . The explicite calculations are given in Appendix where it is stressed that we sho uld be 
distinguish between noninteger and integer S. We apply these results for calculation of eq. (2.21) where 
b{ri) is given by eq.( 2.2C| ).In this case we find 



I{k) ~ wb-— [S{S, z = 1) - Re 3(6, z)] , 



(2.24) 



where S — v{d + 2). Considering non- integer values of the exponent 5 (this assumption will be checked 
afterwards) eq. (|5.2|) yields 



_ d + 2 

I{k) = wb 2 



2p 



(2.25) 



This equation gener alizes eq. ( ^.22| ) ( see also eq. (B.9) of ref. g ). The Euler equation (|]l|) together 
with the expansion (2.25) gives, in the fc limit 



^ 2 ^ 2 2r(<5)cos(f(5) 



0{k% 



(2.26) 



where 5 = v{d + 2). 

Let us impose first the balancing conditions in the spirit of des Cloizeaux: (i) the first two terms on the 
r.h.s. balance each other, (ii) the term on the l.h.s. and the third term on the r.h.s. are also balanced. 
Then we immediately obtain 

I. = 3 , (2.27) 



as a valid exponent, but together with the additional conditions 

wb-^C{5-2) = 1, 
9- 



, _±t2 



2T{5) cos(f^)' 



(2.28) 
(2.29) 



contained in the above equations, the physical interpretation becomes, however, more complicated. More- 
over, we must add to these equations the general relation between g and b, which originates from equation 

o 



1 



2gr(l + 2^) s\n{-Kvy 



(2.30) 



5 



Equation ( ^.27 ) is thus the well known result for the exponent w hereas equations ( 2.28| )-( 2.30 ) set up 
three equations for the two unknowns &, g. Equations ( ^.28 )-( p.30 ) determines not only the amplitudes 



g and h but also fixes the value of the interaction parameter w itself. This delicate point, (see equation 
(IV. 52) of reference [||) has already been discussed by AUegra This is contradictory in the sense 
that generally the interaction parameter w is an arbitrary quantity. The proper balance condition must 
consist of a single equation, whereas the conditions (2.28)-( 2.29p , derived in the asymptotic analysis of 
des Cloizeaux, overconstrains the parameters involved. Besides that the parameter 5 = 2(d + 2)/d is only 
non-integer for d = 3 and is integer for d = 1,2,4 (where the Flory scaling produces exact results). This 
indicates that the des Cloizeaux solution can not be valid in this later case since it relies on the expansion 
( 2.25 ) which (see Appendix) is only valid for non-integer values of i5. 

In order to avoid the overconstraints we try to impose the Flory bala ncing conditions : (i) the elastic 
term, which is proportional to , is balanced by the interaction term eq.(2.21), (ii) the term on the l.h.s. 
(or the entropic term) is negligible compare to the elastic and interaction terms on the r.h.s. Obvio usly 
the Flory balancing conditions have one constraint less. In this case for S{5,z) we must use eq.(5.5) 
where 5 = m ~ vi^d -I- 2) = 3 and z = exp(zfc) . Th is immediately leads to the well known Flory exponent 
vy = 3/(d-|- 2) . Then the interaction term (2.21) becomes 



I{k) = lim wh-— [S{m ^ 3, z = 1) - Re S{m = 3, z)] = 

wb-"^ £ ^ [^,(3) - v.(i) - logfc] ^ + o{k% 



(2.31) 



The resulting Euler equation takes the form 



(2.32) 



Eq.(2.32) can be also obtained directly from eq.(2.26) if we assume 5 = 3-|-e, where e ^ 0. Then the poles 
which have its origin from the second and third terms on the r. h. s. are canceled (see e.g. Appendix 
and sec. 1.11 in ref. ||l^) but the term fc^ logfc appears. The presence of the fc^ logfc term in equation 
( 2.32 ) leads to the conclusion that the pure power law ansa tz g{k) oc k^'^^^ (or h{n) oc in?'^) does not 
solve in the asymptotic limit fc ^ the Euler equation (2.18) once the expon ent 5 = (d + 2)v is assumed 
to be integer. In fact, the direct numerical solution of the Euler equations ( ^.15| ) and ( ^.18| ) presented 
in the next section does not show a pure power law behavior in the asymptotic regime fc << 1. We 
also see that the term fc^ logfc naturally appears in the previous expressions for the series S{5, fc), which 
suggests to consider logarithmic corrections in the initial ansatz for h{n). This will allow to stay within a 
Flory-like solution (corresponding to the integer value 5 = v{d+2) = 3) and also to compute explicitly the 
exponent that characterize the logarithmic correction to the pure power law scaling behavior. Moreover 
the problem of having three equations for the two parameters g and &, which invalidates the asymptotic 
analysis of des Cloizeaux, is naturally solved within the new ansatz we propose in the next subsection. 
The value of w is no longer constrained to a certain value as it should be the case in an ideal experiment 
where it is allowed to tune the Coulombic interaction (or the quality of the solvent for the self avoiding 
problem). 



C. Corrections to scaling 

As it was suggested we will stay with in the Flory balancing condition but try the next asymptotic 
ansatz for the correlation function ( 2.12| ) 

6(n) ~ 6o"^''(logn)^'', (2.33) 

where 7 is the exponent that char acterize the logarithmic correction to scaling. Substituting the ansatz 
(2.33) in the Euler equation ( ^.18 ) gives for the interaction term 

-i+l-^ 1 — cos(nfc) ,„„,x 
I{k)^wbr. ^ y , ^ ^ . 2.34 

The largest contribution to this series at small fc comes from the large n limit. We can substitute the 
sum with an integral (for the power law ansatz, this is explained and justified in the Appendix) to obtain 



6 



m 



dn 



1 — cos(nfc) 



„i.(d+2)[lQg„]7(d+2) ■ 



(2.35) 



where n is a cutoff we introduce to regularize the integral. Equation ( 2.35| ), after substituting n = x/k 
reads 



I{k) ~ ^ 



dx- 



1 — cos(a;) 



(2.36) 



where, because of the Flory balancing conditions (/(fc) oc fc^), we should impose v{d + 2) = 3. Let us 
consider a v alue t such that at x < t it is ensured that 1 — cosx ~ x^/2 holds. In this case the integral 
in equation ( 2.36| ) becomes 



J{k) 



dx 



1 — cos(a;) 



a;3[loga; - log fc]T('^+2) 



dx- 
dx- 



1 



1 

2 Jt,^"^-^ x[\ogx ~\ogk]^id+^) 
1 — cos(a;) 



= Ji + J2 



a;3[loga; - logfc]T('^+2) 
At fc ^ J2(fc) 0, while the integral Ji(fc) defined above can be evaluated as 

1 



Ji{k) 



2[7(d+2)-l] 



log(n)] 



l-7(<i+2)l 



(2.37) 



(2.38) 



provided that 7(d + 2) > 1 and n > 1. We conclude that the Euler equation ( 2.1§| ) is properly regularized 
by the logarithmic ansatz (2.33) and the cutoff parameter n > 1. The balance in equation ( 2.18| ), between 
the elastic energy term and the interaction term, gives 



7 > 



d + 2' 
I 

d + 2' 



(2.39) 
(2.40) 



D. Correction to scaling - charged chains 



For the Coulombic case, considering the dimensional shift rule (2.17), we obtain 



7 > 



3 
1 



(2.41) 
(2.42) 



The inequality we just derived will be verified numerically in the next section. Values for the exponent 



7 will be obtained via the numerical solution of the Euler equations (2.10)-(2.1f:) and via Monte Carlo 
simulations. Within our asymptotic analysis, we can obtain an expression for the exponent 7 that depends 
on the interaction strength. In particular by balancing the prefactors of the proper terms we have 



7 



1 



w 



(2.43) 



where we also took into account that the amplitude 60 can depend from w. The exponent 7 depends 
on the interaction strength, being non univ ersal. This will be tested against simulations and numerical 
minimization of the variational free energy ( 2.1C ). 

It was suggested Q that above a certain minimum charge fraction / ~ iV^'^/'*(6//B)^^^ the electrostatic 
interactions are relevant and the end-to-end distance can be written in the form 



i?~ 7V/2/3(ZB62)(logiV)^. 



(2.44) 
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In this equation the value of 7 = 1/3 is confirmed also by a simplified chain under tension variational 
method |^,^ where the trial Hamiltonian is the Hamiltonian of a Gaussian chain subject to an external 
tension F. Indeed we also reproduce this result for large enough values of the Coulomb strength as 
we will discuss in the next section. We observe that the exponent 7 is non universal, i.e. decreasing 
with the interaction strength, and a pproa ching the limiting value 7=1/3 from above in the region of 
/3 > 1, in agreement with equation ( 2.42 ). At /3 — 1, the value of 7 ~ 0.38 is in good agreement with 
the calculation of Soderberg and collaborators It is important to remark that if we assume in the 
previous asymptotic analysis that the 7 exponent is universal a spurious nested logarithmic form for the 
end-to-end distance (the way AUegra Q discussed long time ago) is obtained. This scenario however is 
not confirmed by the numerical analysis of the Euler equations as will be presented in section III. The 
balancing between the elastic energy and the interaction term determines the exponents (2.39)-( 2.4ClD for 
the self avoiding chain (or shor t-range interaction) problem, whereas for the Coulombic case the result is 

given in equations (2^)- ( 2.42 ). 

The term that appears on the l.h.s. of equation ( |2.26D is negligible in this case. The proper condition 
for this to happen reads \ + 2v > 2, which leads to the condition d < 4 for the self-avoiding chain ||l| 
and d < 6 for the Coulombic chain |2^]. Accordingly the upper critical dimensionality is dc = 4 for the 
former case and dc = 6 for the latter one. At d > dc we can neglect the interaction term and the chain 
becomes Gaussian. 

We would like to make a general remark at this point. It is clear that the variational technique cannot 
cover the results of the RG. From the universality of the self avoiding walk problem it is well known 
that the choice of the parameter w does not change the scaling exponent and below the upper critical 
dimension no logarithmic corrections are necessary. At the critical dimension, these become obviously 
real. On the other hand the variational technique produce these corrections naturally and works well 
in the mean field limit. For the case of polyelectrolytes the situation is more complex. The Coulomb 
potential changes its nature from long range to short range between three and four dimensions. Therefore 
the fixed point and scaling behavior is completely different, as it has been discussed already by de Gennes 
M. It is therefore also necessary to check our results by alternative methods. 



III. NUMERICAL SOLUTION OF THE EULER EQUATIONS 



In this section we consider in detail the Euler equations ( ^.11 ) - ( 2.15| ) for a polymer chain of length 



N with cyclic boundary conditions, obtained by minimization of the variational free energy Fy- Two 
different ways in order to solve the variational principle minimization for a finite chain of length N are 
possible. On the one side, a multidimensional minimization of the free energy Fy in a N dimensional 
space is considered On the other side, a direct solution of the Euler equations will be presented. 

We checked first that our multidimensional algorithm is able to reproduce the simple Gaussian profile 
for the simple case of the non interacting chain problem. We then implement our algorithm by choosing 
a very small value of the interaction strength (3, where the solution is expected to be very close to the 
simple Gaussian chain solution, and give iV-|-l trial profiles for the function gN{k). The multidimensional 
algorithm via a sequence of reflections and translations between the TV -|- 1 trial profiles (k) , relaxes 
to the true solution that can be resolved within the numerical precision desired. Once the algorithm was 
properly tuned, we moved on studying systematically the shape of the optimal profile that minimizes 
the free energy Fy for increasing values of the interaction strength For any specifi c valu e of f3 the 



optimal profile was checked to satisfy the basic symmetry requirements given in equation (2.14). Once the 
numerical solution gN{k) that minimize the variational free energy was found, we proceeded to analyze 
the properties of the solution in the asymptotic regime, in order to check our results for the maximal mean 
square distance b{N/2). Results for the mean square distance {{n+n — ^if) for intermediate values of n 
were also studied in order to confirm the p henom e nolog ical predictions of the electrostatic blob picture. 



A direct implementation of the equations (2.11)- (2.15) does confirm the multidimensional minimization 



of the free energy Fy and also represents a convenient method to solve the variational problem itself. 
We also consider, for the case of A = 1 a Fast Fourier Algorithm that let us consider chains of length 
up to = 2^^. The difference between the values of the exponent 7 computed with N < 1024 and 
N < 32768 are very small and tends to vanish when the interaction strength increases. We consider this 
as a second, independent test to confirm the non-universal behavior of 7. For this reason we will consider 
shorter chains only for what concerns the case of A = 2 and for the Self Avoiding Chain Problem (see 
s ectio n lie) . We start with a purely Gaussian profile in order to compute the left hand side of equation 
( ^.15 ), via the Fourier transform of (2.13) that defines the mean square distance b{n) = {[vi+n — t^Y) 



for a small value of the interaction strength (3 and iterate this procedure until convergence is achieved. 
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Once the solution of the Euler equations is found for a specific value of the interaction, one can use it as 
a starting point for the following iteration procedure, in order to solve the same set of equations for a 
higher value of (3. It is important to note that a simple recursive implementation of the Euler equations 
above is highly unstable in the sense that the number of iterations required as a function of (3 rapidly 
increases, making a simple iteration procedure quite ineffective. A simple solution to this problem can be 
found borrowing ideas from the renormalization group literature . Renormalization Group Theory was 
successfully applied to determine the phase diagram and magnetic properties of disordered spin systems 
. A functional set of recursion relations for the probability distribution function of bond randomness 
was studied and a fixed distribution was found for different values of the initial conditions. Instead 
of a simple iterative procedure, which could be a rather slow converging process, a multidimensional 
bisection between different initial conditions for the renormalization group trajectory turns out to be a 
very effective method in order to achieve convergence. In a similar fashion, instead of a simple iteration 
of the Euler equations above, we prefer to bisect between possible initial profiles of the function g{k) 
given above. 



A. Results for the Coulombic chain 



We obtained the solution of the Euler equations ( 2.11 )-( 2.15 ) via the iterative method discussed above 
and independently via the direct minimization of the variational free energy Fy in equation (2.1C) within 
the numerical precision . We considered cyclic chains of length 2^ < < 2^^ — 32768, for different 
values of the interaction strength (3 and two different forms of the potential (2.8), corresponding to the 
values A = 1,2. We report the results of the end-to-end distance for different values of A and for chains 
of increasing size A^ at values of /? = 0.5 and (3 — 1.0. 



A 


/? 


N=120 


N=240 


N=320 


N=520 


N=1024 


N=2048 


N=4096 


N=8192 


N=16384 


N=32768 


1.0 


0.5 


60.6 


128.2 


174.5 


293.2 


601.7 


1249.4 


2583.4 


5372.2 


11071.2 


22769.8 


1.0 


1.0 


78.0 


164.4 


223.6 


374.2 


767.59 


1590.6 


3283.9 


6778.8 


13942.4 


28627.3 



A 


/3 


N=120 


N^240 


N==320 


N==520 


N=1024 


2.0 


0.5 


26.2 


48.6 


63.0 


98.12 


140.2 


2.0 


1.0 


33.2 


62.5 


81.4 


126.2 


248.7 



Tah.l: Results for the end-to-end distance R for different values of the chain length A^ an d the interaction 
strength (3 and A. These results are obtained solving explicitly the the Euler equations (2.11)-(2.15). 

Let us begin to discuss the case of a simple polyelectrolyte chain in d = 3 with pure Coulombic 
interaction (A = 1). We analyzed the end-to-end distance behavior for increasing values of A^, i.e. 
h{N/2). This is reported in Fig.l. Dependence of the exponent 7 on the interaction strength suggests a 
weak violation of universality as manifest in Fig. 2. 
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Chain Length, N 



FIG. 1. End-to-End distance R versus the chain length A*', for increasing values of the interaction strength 
/3. The inset shows the logarithmic plot, when a linear dep endence of R on the chain length is assumed, beside 
logarithmic corrections, for values of A = 1 in equation (2.8). 
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FIG. 2. The exponent 7 as measured (see Fig.l) for different values of the interaction strength /3. The exponent 
7 is seen to decrease for large values of the interaction strength, towards the expected value 7 — 1/3. The squares 
indicates the computed value of 7 with chains of length A'^ < 1024 while stars show the computed value of 7 for 
chains of length A'' < 2^^. 



The des Cloizeaux argument, once extended to the polyelectrolyte chain problem in d = 3 imphes an 
overstretching {u = 2/A = 2), while the ansatz suggested above indicates v = 3/(A + 2) = 1. In order 
to compute the non- universal exponent 7 we plot J?/jV '^/(-'^+^) as shown in the inset of Fig.l. The non- 
universal behavior of 7 agrees with the inequahty ( 2.42| ) derived in the asymptotic analysis of section II. 
In Fig. 2, the exponent 7 decreases with increasing values of /3, approaching values very close to 7 = 1/3, 
for large values of the interaction strength. A recent discussion about the logarithmic correction to scaling 
behavior for polyelectrolyte chains, has been given p9| and a direct derivation of the value 7 = 1/3 has 
been obtained. In the notation of reference psj ] the small temperature limit T — > corresponds to the 
large interaction strength limit /3 — > 00 which is consistent with our findings. The result 



R(x N{\ogNy/^ 



(3.1) 
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obtained in the large (3 ^ oo limit, can be interpreted as a stretching of the harmonic bonds of the chain. 
The value of 7 = 1/3 was also obtained by a simplified chain under tension variational method where 
the trial Hamiltonian is the Hamiltonian of a Gaussian chain subject to an external tension F. This 
result, confirmed also by the Monte Carlo simulations presented in the next section, holds in the limit of 
large interaction strength but should not considered as a limiting case of little practical impact. Ordinary 
temperatures in the aqueous solutions correspond to surprisingly "high" values of the interaction strength 

A rather important problem discussed in the past H is the polyelectrolyte chain with an interaction 



(2.8) at A = 2. After the proper rescaling of the interaction strength /3 (see equation 2.10), this problem 
can be thought as a chain with pure Coulombic interaction for d = 4 . Results for the end-to-end distance 
versus the chain length iV are shown in Fig. 3. 
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FIG. 3. End-to-End distance R versus the chain length TV for increasing values of the interaction strength /3 
at A = 2 . The inset shows the logarithmic plot, when a power law A'^'^^ dependence of R on the chain length is 
assumed, beside logarithmic corrections. 

A similar discussion about the non-universal behavior of the exponent 7, follows the same lines as 
in the previous case (Fig.l-Fig.2). According to the ansatz we proposed in the asymptotic analysis of 
section II C we compute the exponent 7 as in Fig. 4, where — vp is the Flory exponent. 
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FIG. 4. The exponent 7 as measured from Fig. 3 for different values of the interaction strength /3. The exponent 
7 is seen to increase for large values of the interaction strength towards the value 7 ~ 1.0. 
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In Fig. 5 we compute the non-universal behavior of the exponent 7 assuming a linear dependence of 
the end-to-end distance of the chain length TV, beside logarithmic terms. This is not the case, for the 
potential (|]^) with A = 2, but interesting comparison with previous results follows. The exponent 7, 
as in Fig. 5, approaches values of 7 « —0.4 at values of /3 = 1.0 where Monte Carlo calculations were 
performed |23] and where the computed value of 7 « —0.35 is in good agreement with our result. A two 
parameter fit would represent a very conclusive way to discriminate between these two possibilities, but 
the chain lengths that one can reach in simulations and in a direct minimization of the variational free 
energy are not long enough to make such analysis effective. 
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FIG. 5. The exponent 7 as measured from Fig. 3 for different values of the interaction strength /3. The exponent 
7 is seen to increase for large values of the interaction strength towards the value 7 « 1.0, once a linear dependence 
of R on the chain length is assumed, beside logarithmic corrections, for values of A = 2 in equation (2.8). This 
results should be compared with the value of 7 ~ —0.35 previously found via Monte Carlo simulations, for values 
of /? = 1.0. 



B. Electrostatic blob picture 



We will now present a numerical analysis of the electrostatic blob picture. We measured the crossover 
(shown by the mean square distance as a function of the polymer chain length N) from a Gaussian 
regime, observed at small distances, to the stretched regime for long enough distances. This represents 
an effective way to compute the electrostatic blob size and to check our method. One of the major 
achievement of the Gaussian variational approach discussed in section I is that it does not simply provides 
a quantitative description of th e large n behavior correlation function. A full numerical solution of 
the Euler equations ( ^.11 )-( ^.15 ) provides informations on the intermediate behavior of the correlation 
functions 6(n) = ((ri+„ — t^)'^). That is why it is interesting to check whether short enough segments of the 
polymer chains behave as Gaussian chai ns, as expected from general scaling arguments whereas the 
large n behavior corresponds to equation (2.44). Once the crossover is observed, the size of an electrostatic 
blob ^ can be extrapolated. The dependence of the electrostatic blob size ^ from the interaction strength 
/?, that we expect from phenomenological arguments, is indeed confirmed within our variational approach 
(see Fig. 6). In particular we expect 



(3.2) 



In Fig. 6 we show the electrostatic blob size as a function of the interaction strength /?. The ^ cx /3 ^/■^ 
behavior expected from standard scaling arguments is here also found. 



12 



ai" - 

N 

_o 
CO 



Fit Parameters A=1.89 
B=1.34 



0.01 0.02 0.03 

Interaction Strength, p 



FIG. 6. Blob size, as measured from the direct solution of equations ( p.l5| )-(2.10), considering the crossover 
between the small distance behavior of the mean square distance b{n) versus the large distance behavior. Different 
values of the electrostatic blob size for increasing values of the interaction strength (3 are compared with the 
f ~ f3~^^^ behavior expected from phenomenological arguments. 



C. Results for the Self Avoiding Chain problem 



We consider the problem of a polymer chain with a short range interaction between the monomers, in a 
good solvent. Within the Gaussia n Va riation al Pr inciple, the expression for the free energy and the Euler 
equat ion are given by equations ( 2.16| ) and ( [2.18 ). In section II we have discussed the dimensional shift 
( 2.17 ). Nevertheless this equivalence is purely formal and no deep insight is hidden behind it. The two 
problems considered are indeed pretty different and they reduce to the same type of variational equations 
because of the Gaussian approximation itself. It is important to remark that the Gaussian Variational 
Principle, specially reliable for the long ranged polyelectrolyte problem, has been originally introduced to 
study the self avoiding homopolymer chain problem |^,^ . Even though it is nowadays well established that 
the Gaussian Variational Principle is not very well suited for the short ranged interacting chain problem, 
we wish to make clear that within this approximation, apart for spurious logarithmic corrections, the 
usual Flory-like power law behavior should be recovered. 

In the original asymptotic analysis of des Gloizeaux the non integer value of 5 is tacitly assumed, so 
that the result oiv = 2/d for the swelling exponent was obtained. Although we believe this result to be 
incorrect, it is often mentioned in the literature to discuss the asymptotic properties of polyelectrolyte 
chains According to our analysis and to the previous analysis of Allegra and coworkers the 

problem of a self avoiding chain in good solvent turns out, within the Gaussian Variational Principle, to 
present the expected Flory-like behavior. This is not surprising since the original derivation of Flory also 
belong to the class of Gaussian Variational approximations, since a Gaussian probability distribution of 
the distance between the ends of the chain is assumed. A numerical justification of the des Gloizeaux 
result for the swelling exponent can also be found in the literature We believe that this numerical 
analysis is inconsistent since the constraint on the interaction parameter (assumed by des Gloizeaux and 
critically reviewed in section II) is enforced in the numerical algorithm. This makes the results obtained 
highly questionable. In Fig. 7 and Fig. 8 we present the results of the end-to-end distance for the self 
avoiding homopolymer chain problem studied via the numerical solution of the Euler equations discussed 
above. 
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FIG. 7. End-to-End distance R versus the chain length N. for increasing values of the interaction strength 
w, that represent in this case the quality of the solvent. The inset shows the logarithmic plot, when a iV"^ 
dependence of R on N is assumed, beside logarithmic corrections. 

Results for the cnd-to-cnd distance as a function of the polymer chain length N arc shown and a 
quantitative analysis of the non universal logarithmic part to be considered as a spurious effect due to 
the short range nature of the problem is presented. 
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FIG. 8. The exponent 7 as measured from Fig. 7 for different values of the interaction strength w. The exponent 
7 is seen to increase for large values of the interaction strength. 

According to the asymptotic analysis given above we expect the exponent 7 to be non universal and 
to depend indeed on the quality of the solvent. 



IV. NUMERICAL SIMULATIONS 



In order to test the conclusions obtained in the previotis section we decided to perform a numerical 
simulation for a polyelectrolyte chain in an ideal salt free solution, when pure Coulombic interactions are 
considered between the monomers. A bit of caution is required since we will compare the results of our 
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simulations with the results obtained via the Gaussian variational principle, where a cyclic geometry for 
the polymer chain is chosen. It was easy however to device a numerical simulation of the Monte Carlo 
type where close chain configurations can be taken under investigation. As remarked long time ago, |^ 
the cyclic geometry chosen within the variational calculation becomes unessential as soon as long enough 
chains are considered and this turns out to be confirmed by our numerical simulations. We then consider 
a Monte Carlo simulation, performed in the canonical ensemble via the standard Metropolis algorithm. 
We consider a pivot move for what concerns the elementary Monte Carlo move, i.e. we assume rigid bonds 
between the monomers and an elementary pivot move consists in a random rigid rotation of a randomly 
chosen part of the chain. The pivot algorithm, first described in reference and lately applied in the 
context of polyelectrolyte chains with and without salt jl^ is a very effective method to confirm our 
general conclusions. Several million Monte Carlo moves (the number which was properly increased with 
the chain length) are required in order to reduce below 1% percent the statistical fiuctuations for the 
end-to-end distance. We also combine our Monte Carlo simulation with a Molecular dynamics subroutine 
that is randomly called, during the Monte Carlo simulation, with a certain fractional occurrence ratio. 
On the one side, the pivot algorithm is well suited to simulate the polymer chain at long enough scales. 
The pivot moves are in a sense collective moves and a good description of the chain is certainly obtained 
at long scales. On the other side the Molecular Dynamic algorithm is well suited to simulate the chain at 
small distances. We initially perform independent simulations with the Monte Carlo and the Molecular 
dynamics algorithm in order to tune the simulation parameters involved for a chain of a given length N. 
Once the two methods were independently tested, we proceeded to combine them as discussed above. 

We now present our results for the end-to-end distance for d = 3 polyelectrolyte problem. In Fig. 9 
the end-to-end distance is shown as a function of the polymer length N for the specific value of P — 1.0. 
The inset shows the end-to-end distance divided by the pure power law Flory-like dependence N'^ as a 
function of log A^. We repeated the same analysis of Fig. 9 changing the value of the interaction strength 
p. The non universal behavior obtained in the asymptotic analysis of section II and in the numerical 
analysis of the Euler equations (2.11) - ( p. 15 ), is found within our simulation results (see Fig. 10). 
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FIG. 9. End-to-End distance R versus the chain length A'^ for the polyelectrolyte chain at increasing values of 
the interaction strength /3 and d — 3, obtained via the pivot Monte Carlo algorithm combined with Molecular 
Dynamics for the value of the occurrence factor r = 0.1. Different simulations, corresponding to different values 
of the occurrence factor r have been shown to be numerically consistent within each other. The inset shows the 
logarithmic plot, when a linear dependence of i? on iV is assumed, beside logarithmic correction 



The non universal exponent 7 is seen to decrease, as we obtained with the variational method of section 
II-III. Nevertheless one should treat the results of simulation with a bit caution. 



15 







A=1.0 
V=1.0 


1 

y 






0.5 







0.4 0.8 1.2 

Interaction Strength, P 



FIG. 10. The non universal exponent 7, for different values of the interaction strength (5 as measured from 
Fig. 9, for the specific value oi (5 = 1.0 and for the values (5 = 0.5, 0.75, 1.5. 

In the numerical simulations monomers are considered to be rigidly connected while in the variational 
principle approach monomers are connected by springs. In this case one can not expect that the exponent 
7 as a function of interaction strength has the same large /? behavior. 



V. CONCLUSIONS 



In the present paper we have reinvestigated the Gaussian variational principle and its applicability to 
polymer chain with long and short range interactions. The asymptotic analysis originally presented by 
des Cloizeaux in the context of self avoiding walks has been critically reconsidered. In particularly we 
have shown that the des Cloizeaux's analysis overconstains the model's parameters. It has been found 
that within the Flory type of balancing conditions (whe n the elastic term is balanced by t he in teraction 
but the entropic term is not relevant) the asymptotic law ( p. 33 ) satisfies the Euler equation ( ^.15 ) . In this 
asymptotic law the Flory exponent v — i / {d -\- 2) and logarithmic corrections to scaling is characterized 
by a non - universal exponent 7. Within the Gaussian variational principle the universal expo nents for 
the case of the Coulombic and short range interactions are related by a dimensional shift ( 2.17 ). 

Gaussian variational principle suits much better for the long range problems which was already dis- 
cussed in reference [fzi. The arguments we have presented prove that the scaling form ( 2.33 ) is val id for 

The exponents v and 7 satisfy the equations (2.41) 
universal and depends from the 

interaction parameter as well as from the chain model (compare e.g. Fig. 2 and Fig. 10) . In this sense we 
may consider the combination l^s = a{\ogN)'^ as an effective Kuhn segment for chains with unscreened 
type of interactions. 

On the contrary, for the short range interactions the Gaussian variational principle simply overestimates 
the probability of the segment - to - segment contacts which leads to more extended configurations 
compare to the Monte Carlo results |2^] . We conclude observing that the models which are characterized 
by the interplay between long - ranged attraction and repulsion are also useful systems to be treated 
by the Gaussian variational principle, as presented here. Appropriate systems are polyampholyte | p6[ |. 

have found strongly nonuniversal 



a gener al long - range interaction pote ntial (2. 
- (2.42). According to the condition ( ^.42 ), the exponent 7 is non 
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For the randomly charged polyampholytes Kantor and Kardar 
behavior. For example at d < 4 the spatial configurations depend on the polyampholyte excess charge Q 
so that if Q > Q* « e^/N the compact form changes to the stretched one. We believe that the Gaussian 
variational principle is a good starting point for the investigation of such systems. 



APPENDIX A 



The asymptotic analysis of equations (2.18) suggests to consider the following series 
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n=l 



(5.1) 



In the asymptotic limit fc ^ the interaction term of equation (2.18) is si mply related to S[5^ z). Even 
the simple power law ansatz h{n) cx n"^ requires some care. The expression ( 5.1) is evaluated exactly but 
two rather different expressions are obtained if the exponent 5 = i'{d + 2) is assumed to be integer or 
non-integer. The case of the exponent S being non-integer reduces the interaction term to the expression 
given in equation ( 2.25 ) and leads to the des Cloizeaux result. Assuming an integer value for the exponent 
S produces instead a different expression for the interaction term and the proper balancing condition with 
the elastic energy naturally follows. 



Let us consider the series S{S,z) in this two distinct cases. In equation (5.1) z and S can be complex 
quantities, as soon as \z\ < 1 and Re 5 > 1 at |z| = 1 whereas i?e (5 > if |z| < 1. The calculation of the 
series (K^) ( see sec. 1.11 in ref. [H ) gives 



SiS,z) 



r(i-<5)[iog(V^ + Ec('5-0-^'°''^' 

z ^-^ 

r=0 



(5.2) 



where |log^| < 27r, T{x) is the gamma function, S ^ 1,2,.. and C,(x) is the C-zeta Riemann's function. 
We ar e interested in the case of z = exp(ifc) at fc ^ 0. Con sider ing the real part of both sides of equation 
( 2.23 ) in the asymptotic limit fc — > 0, we obtain equation ( ^.25|) . 

We now consider the case of integer values of b. The ^{x) has poles at all negative integer arguments 
whereas the pole of C,{x) is placed at a; = 1. We write ^ = rn + e where m is a positive integer and e — > 0. 
Then in the vicinity of poles the gamma and C - functions can be rewritten as 



r(l-m-e) 



f-D™ f 1 



(to 



1 



C(l + e) = ^ - - + 0(e) 



We should also take into account that 

log ' ^ 



= 1 + e log log 



(5.3) 



(5.4) 



After that the equation (5.2), due to the cancellation of poles in gamma and C - functions at small values 
of e, becomes M 

floKz)™"^ 

5(m, z) = V 1M (V'M - V^l) - loglog(l/z)) 



(to - 1)! 

OO 

(logzy- 



(5.5) 



where the prime indicates that the term r = to — 1 is to be omitted. In equations ( p.5| ) ( |5.3| ) ^^ix) is the 
digamma function, defined as the logarithm ic d erivative of the gamma funct ion "^ {S) = in the 



dS 



limit A: — > 0, for z — exp(ifc) and S — 3, eq. (5.5) have been considered in eq.(2.31). 

We want to show now that the integral representation of the sum is correct in the asymptotic limit 
and that finite n < I terms do not contribute in this limit. Let us consider explicitly the presence of a 
finite cutoff in the series (5.1). 



n—l 



(5.6) 



which reduces to equation (5^) at / = 1. For integer values of S it is easy to derive exactly the following 
expression [ p6[ : 

5(m, z, I) = Wm) - m - loglog(l/z)) 



(m - 1)! 

OO 

+ E 



(5.7) 



r=0 
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Consider the series (q^) with a finite cutoff I in the asymptotic Hmit fc — > 0. 
contributions to the series from finite n < I terms, when m = 3. We obtain 



E 

n=l 



1 — cos(nA:) k'^ 



(V;(3)-V(0-logfc)+O(fc4), 



We want to evaluate the 



(5.8) 



so that the contribution from finite n takes the following form 



' 1 



cos(nfc) 



(5.9) 



Summing up equations (5^) and ( |5.9| ) we obtain again the result ( 2.31 ). The expression ( |5.8| ) shows that 
for small k values only large n terms contribute to the sum (p]q). Indeed the contribution coming from 
n > I produces a term log(fc). This, whenever k << exp(— ?/;(?)), will dominate the simple behavior. 
Moreover, for large values of I (see ref. [^) we have ip{l) oc log(Z) and the condition above takes the form 
k << l/l. It follows that for k << l/l the term fc^ log(fc) dominates and the sum does not depends on 
the low limit value. In this case, the series can be well approximated by an integral, i.e. 



^(fc) = E 

n—l 



1 



i{nk) 



dn 



1 — cos(nfc) 



= k' 



dx 



1 — cos(a:;) 



= k^I{k) 



(5.10) 



kl 



In order to compute I{k), consider the derivative of the previous expression, I'{kl) = ~l/2kl so that 
I{kl) ~ —J log(fcZ) + c, where c is a constant, and S{k) = ^{c — log{kl)) + 0{k'^). Comparison with 
equation (^^) shows that for k « l/l only large n contributions are relevant. 
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